package simuspec;
import peptide.Peptide;
import spectrum.Spectrum;
public class SimuSpect {
	Peptide pep;
	Read_Weca_Result wr;
	int mer_size=2;
	int Yup = 7;
	int Bup=7;
	public static double[][]peakGroup;
	double initial_Intensity=10000;
	public SimuSpect(String pepStr,String wr){
		this.pep = new Peptide(pepStr,pepStr.length());
		//if(this.pep.)
		this.wr = new Read_Weca_Result(wr);
	}
	
	public Spectrum get_theor_spectrum_with_iso(){
		//System.out.println("get_theor_spectrum_with_iso()");
		
		double[][] intensity_array_y =this.get_theor_yArray_with_iso();
		//publicMethod.ArrayManipulation.printArray(intensity_array_y);
		double[][] intensity_array_b =this.get_theor_bArray_with_iso();
		//publicMethod.ArrayManipulation.printArray(intensity_array_b);
		double     precusor_mass=this.pep.getPrecursorMass();
		int  charge=2;
		double fragmentError=0.1;
		double precusorError=3.0;
		double intensityThres=0;
		double[][] int_a =publicMethod.ArrayManipulation.array_merge(intensity_array_y,intensity_array_b);
		//publicMethod.ArrayManipulation.printArray(int_a);
		int_a = publicMethod.ArrayManipulation.ascendArray(int_a);
		Spectrum spe = new Spectrum(int_a,precusor_mass, charge,fragmentError, precusorError,intensityThres);
		//spe.print();
		return spe;
	}

	
	public Spectrum get_theor_spectrum_y_with_iso(){
		
		double[][] intensity_array_y =this.get_theor_yArray_with_iso();
		double     precusor_mass=0;
		int  charge=0;
		double fragmentError=0.1;
		double precusorError=3.0;
		double intensityThres=0;
		double[][] int_a = publicMethod.ArrayManipulation.ascendArray(intensity_array_y);
		Spectrum spe = new Spectrum(int_a,precusor_mass, charge,fragmentError, precusorError,intensityThres);
		//spe.print();
		return spe;
	}
	
	public Spectrum get_theor_spectrum_b_with_iso(){
		
		double[][] intensity_array_b =this.get_theor_bArray_with_iso();
		double     precusor_mass=0;
		int  charge=0;
		double fragmentError=0.1;
		double precusorError=3.0;
		double intensityThres=0;
		double[][] int_a = publicMethod.ArrayManipulation.ascendArray(intensity_array_b);
		Spectrum spe = new Spectrum(int_a,precusor_mass, charge,fragmentError, precusorError,intensityThres);
		//spe.print();
		return spe;
	}



	public double[][] get_theor_bArray_with_iso(){
		
		double[][] intensity_array_b =new double[pep.pep_length*3-6][2];//peptide intensity array;k 
		String skel= pep.skeleton_string;
		double inten=this.initial_Intensity;
		for(int i=2;i<skel.length();i++){
			double predict_ln_ratio_YY=0;
			double predict_ln_ratio_BY=0;
			int index = i;
			if(i>=Bup){
    			index=Bup;
 			}
			predict_ln_ratio_BY= get_ln_ratio(i,wr.adjust_factor_BY[index],wr.factors_AA_BY[index],wr.factor_b_BY[index],
					wr.factor_Nterm_BY[index],wr.factor_K_Cterm_BY[index],wr.factor_R_Cterm_BY[index]);
			double inten_b=inten*Math.exp(predict_ln_ratio_BY);
			intensity_array_b[(i-2)*3][0]=pep.getBMass(i);
			intensity_array_b[(i-2)*3][1]= inten_b;
			double[] iso_B = pep.getB_ion(i, 1).get_isotope_ratio();
			intensity_array_b[(i-2)*3+1][0]= pep.getBMass(i)+1;
			intensity_array_b[(i-2)*3+1][1]= inten_b*iso_B[1]/iso_B[0];
			intensity_array_b[(i-2)*3+2][0]= pep.getBMass(i)+2;
			intensity_array_b[(i-2)*3+2][1]= inten_b*iso_B[2]/iso_B[0];
			index = i;
			
			if(i>=Yup){
    			index=Yup;
    			
			}
			predict_ln_ratio_YY = get_ln_ratio(i,wr.adjust_factor_YY[index],wr.factors_AA_YY[index],wr.factor_b_YY[index],
					wr.factor_Nterm_YY[index],wr.factor_K_Cterm_YY[index],wr.factor_R_Cterm_YY[index]);

			inten = inten/Math.exp(predict_ln_ratio_YY);
			
		}
		double[][] int_a =publicMethod.ArrayManipulation.ascendArray(intensity_array_b);
		return int_a;
	}
	public double[][] get_theor_yArray_with_iso(){
		
		double[][] intensity_array_y =new double[pep.pep_length*3-3][2];//peptide intensity array;k
		String skel= pep.skeleton_string;
		double inten=this.initial_Intensity;
		for(int i=2;i<skel.length();i++){
			double predict_ln_ratio_YY=0;
			int index = i;
			
			if(i>=Yup){
    			index=Yup;
    			
			}
			predict_ln_ratio_YY = get_ln_ratio(i,wr.adjust_factor_YY[index],wr.factors_AA_YY[index],wr.factor_b_YY[index],
					wr.factor_Nterm_YY[index],wr.factor_K_Cterm_YY[index],wr.factor_R_Cterm_YY[index]);
			intensity_array_y[(i-1)*3][0]= pep.getYMass(i);
			intensity_array_y[(i-1)*3][1]= inten;
			double[] iso_Y = pep.getY_ion(i, 1).get_isotope_ratio();
			
			intensity_array_y[(i-1)*3+1][0]= pep.getYMass(i)+1;
			intensity_array_y[(i-1)*3+1][1]= inten*iso_Y[1]/iso_Y[0];
			intensity_array_y[(i-1)*3+2][0]= pep.getYMass(i)+2;
			intensity_array_y[(i-1)*3+2][1]= inten*iso_Y[2]/iso_Y[0];
			
			inten = inten/Math.exp(predict_ln_ratio_YY);
			
		}
		double[][] int_a =publicMethod.ArrayManipulation.ascendArray(intensity_array_y);
		//spe.print();
		return int_a;
	}
	
	public double get_ln_ratio(int pos,double[] factor, double[][]  factors_AA,double factor_b,double factor_Nterm,double[] c_K,double[] c_R){
		double ratio=factor_b;
		//publicMethod.ArrayManipulation.printArray(factors_AA);
		char[] AA_array = pep.skeleton_string.toCharArray();
		//AA_array = chaos_array(AA_array);
		for(int i=0;i<AA_array.length-1;i++){
			int y=Character.getNumericValue(AA_array[i]);
			int x=i-pos+Paras.zero_pos;
			if(x<0) x =0;
			//System.out.println(y+"\t"+x);
			ratio += factors_AA[y][x];
			if((AA_array[i]=='K'||AA_array[i]=='R')&&factors_AA[y][x]==0){
				ratio = Double.NaN;
			
			//System.out.println(pep.pep_string+"\t"+pos+"\t"+AA_array[i]+"\t"+x+"\t"+factors_AA[y][x]+"\t"+ratio);
			}
		}
		
		int x=AA_array.length-1-pos+Paras.zero_pos;
		if(AA_array[AA_array.length-1]=='R'){
			ratio += c_R[x];
			//System.out.println(c_R[x]+"\t"+x);
		}
		else if(AA_array[AA_array.length-1]=='K'){
			ratio += c_K[x];
		}
		else{
			int y=Character.getNumericValue(AA_array[AA_array.length-1]);
			ratio += factors_AA[y][x];
			
		}

		double Nterm=0;
		if(pos>mer_size){
			Nterm = pep.getBMass(pos-mer_size);
		}
		ratio += factor_Nterm*Nterm;
		//System.out.println("factors_b\t"+factor_b);
		//System.out.println("factors_N\t"+factor_Nterm+"\t"+Nterm);
		return ratio;
		//return ratio*factor[0]+factor[1];
	}
	public void setPeakGroup()
	{
	   Spectrum thero_sp=get_theor_spectrum_with_iso();
	   peakGroup= thero_sp.intensity_array;
		
	}
	public static void main(String[] args)
	{
		String pepStr="GPAAIQK";
		String weka_file="c:/pa/weca.result";
		
		SimuSpect simu=new SimuSpect(pepStr,weka_file);
		Spectrum thero_sp=simu.get_theor_spectrum_with_iso();
		double[][] peak_group=thero_sp.intensity_array;
		System.out.println(thero_sp.precusor_mass+"\t"+thero_sp.charge);
		for(double[] peak:peak_group)
		{
			System.out.println(peak[0]+"\t"+peak[1]);
		}
		
	    
	}
	
}
